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APPLICATION  OF  HARMONIC  ANALYSIS  METHOD 
TO  RESEARCH  ON  ROTOR  AIRLOADS 

Qiu  Zhenhan 

( Chinese  Helicopter  Research  and  Development  Institute) 

Abstract 

According  to  the  rotor  vortex  theory*  the  rotor  circulation  and  the 
rotor  induced  velocity  are  developed  into  Fourier  series;  The  circulation 
distribution  along  blade  spanwise  is  expressed  in  Terms  of  segment-by¬ 
segment  linear  functions.  In  consequence  the  induced  velocity  equations 
and  the  circulation  equations  are  derived.  The  engineering  application  of 
the  rotor  vortex  theory  is  provided.  Then  the  induced  velocity  and  its 
harmonic  components  are  obtained  to  provide  a  quantitative  basis  for  the 
vortex  model.  For  calculating  each  order  harmonic  components  of  the 
induced  velocity  a  simplified  method  is  put  forward  which  considers  the 
effects  of  each  order  circulation  with  neglecting  those  of  higher  order. 
The  method  saves  the  computer  time  and  is  of  significant  benefit. 


APPLICATION  OF  HARMONIC  ANALYSIS  METHOD  TO  RESEARCH  ON  ROTOR 
AIRLOADS 


Qin  Zhenhan 

Chinese  Helicopter  Research  and  Development  Institute 
Submitted  31  July  1984 


This  paper  uses  the  rotor  vortex  theory 
of  Professor  Wang  Shichuen  as  a  basis  to  develop 
the  rotor  circulation  and  induced  velocity  into 
Fourier  series.  The  circulation  distribution 
along  blade  spanwise  is  expressed  in  terms  of 
segment-by-segment  linear  functions  and  the 
linear  equations  for  calculating  rotor  airloads 
using  the  circulation  as  the  unknown  are  derived. 
This  has  resolved  the  practical  application  prob¬ 
lems  of  the  rotor  vortex  theory  that  for  many 
years  have  not  been  able  to  be  resolved. 


I.  Preface 


There  have  been  many  theories  and  articles  on  the  aerodynamic 
calculations  of  helicopter  rotor-blades,  and  among  which  is  the 
rotor  vortex  theory  of  the  stationary  vortex  system  of  Professor 
Wang  Shicun^1^.  But  it  requires  proper  modifications  when  applying 
this  theory  to  engineering  calculations.  Due  to  the  limitation 
of  resources  in  the  past,  the  rotor-blade  could  only  be  treated 


as  a  rigid  blade,  and  this  is  obviously  improper.  As  the  application 
of  computer  develops  in  our  country,  the  elastic  deformation  (bending, 
twisting,  etc.)  of  the  rotor-blade  can  be  considered  simultaneously 
and  this  can  more  truthfully  reflect  the  rotor  characteristics. 

According  to  Reference  [1],  this  paper  develops  the  rotor  circula¬ 
tion  and  induced  velocity  into  the  Fourier  series.  This  manipulation 
makes  the  physical  significance  between  the  circulation  and  its 
related  parameters  become  very  clear  and  is  rather  conducive  toward 
the  understanding  of  the  vortex  model  and  the  effects  of  various 
vortex  systems.  Although  the  mathematical  manipulations  are  quite 
complex,  a  smaller  computer  capacity  is  required.  While  high  accuracy 
computations  can  be  conducted  on  a  computer  with  large  capacity, 
results  can  also  be  obtained  using  a  computer  with  small  capacity. 

For  example,  it  only  requires  half  an  hour  to  obtain  results  that 
meet  engineering  requirements  on  the  DJS-21  computer  which  only 
processes  ten-thousand  operations  per  second.  Therefore,  it  is 
economical  and  suitable  for  engineering  applications. 

II.  Induced  Velocity  Equations 

According  to  Reference  [1],  the  axial  induced  velocity  at  any 

point  on  the  rotor  hub  can  be  developed  into  the  Fourier  series 
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Here  every  induced  velocity  is  agitated  by  the  circulation  which 
consists  of  attached  vortex,  vertical  free  vortex  and  horizontal 


free  vortex. 


According  to  Reference  [2],  if  the  circulation  distribution 
along  spanwise  is  expressed  in  terms  of  segment-by-segment  linear 
functions,  the  induced  velocity  derived  from  vertical  free  vortex 


can  be  expressed  as 
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and  Cj[m]  are  the  quantites  related  to  Kv  and  the  harmonic  order 
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C.  is  the  quantity  related  to  the  integration  of  the  super  geometrical 
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function 
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For  example 

m,  («  +  1  )/2,  ( m  +  l  )/2,  F,  0) 

The  expressions  of  induced  velocities  v_  and  v„  agitated  by  the  attach- 
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ed  vortex  and  horizontal  free  vortex  are  similar  to  Equation  (2). 


III.  Circulation  Equations 


According  to  Reference  [2],  the  circulation  is  developed  into 
Fourier  series 
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The  circulation  distribution  along  spanwise  is  expressed  in  terms 

of  segment-by-segment  linear  functions,  then  the  circulation  equation 

in  the  form  of  T  is 
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B  is  the  quantity  related  to  B  ,  B  and  motion  parame  >rs.  cor 
ns  n  2  o  s  ^ 

example  5„— n(5,+5,)/{  2  C(5,/3)(  1  —  ?*)+(^V2)(  1  —  ?,)5,)' 
where  5,-l/C  l  +AJCr')i  5,-1/Cl  +AJ Cr(  1  +C*)D 

Flop'  Pisp  and  ot^er  circulation  equations  can  be  derived  according¬ 


ly.  Rotor-Blade  Airloads 


Each  order  harmonic  component  of  circulation  can  be  obtained 
from  the  linear  algebraic  equations  with  circulation  as  the  unknown 
which  are  derived  by  substituting  the  induced  velocity  equations  into 
the  circulation  equations.  According  to  the  Zhukovskiy  formula,  the 
airloads  per  unit  rotor-blade  length  are: 
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Thus,  the  expressions  of  each  order  airloads  per  unit  rotor-blade 

length  are  obtained;  and  therefore  the  airloads  are  obtained. 
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V.  Examples  and  Discussions 


Numerical  computations  were  conducted  using  the  H-34  helicopter 
as  an  example.  The  calculated  values  were  compared  with  the  measured 
values  during  actual  flight  and  the  results  are  satisfactory;  see 
Figs.  1  and  2. 


Another  helicopter  with  takeoff  weight  of  14,400  kg  and  rotor 
diameter  f  21.3  m  was  used  as  an  example  to  calculate  each  component 
of  the  induced  velocity  and  obtain  the  curves  of  the  azimuthal  varia¬ 
tion  of  total  axial  induced  velocity '( Fig .  3).  Due  to  the  limitation 
in  space,  only  the  second  order  components  of  the  induced  velocities 
which  are  agitated  by  separate  harmonic  components  of  attached  vortex, 
vertical  free  vortex  and  horizontal  free  vortex,  are  given.  See 
Table  1-3.  The  obtained  harmonic  components  of  the  induced  velocities 
have  not  yet  been  seen  in  any  references.  They  provide  a  quantitative 
basis  for  the  analysis  and  research  of  the  vortex  model. 
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Comparison  of  the  azimuthal  variation  of  “77“*(N/m)  at  U  "0.2,  for  * » 0 

- terti  — - calculation . 


rig.  3  .  Azimuthal  variatioa  of  the  induced  velocity  v  at  I*  ■  0 . 30.  for  *  “  0 . 75 

I  —The  induced  velocity  c/»  by  bound  vortices.  I  — The  induced  velocity  pj  by  longitudinal  vortices. 

I— The  induced  velocity  v»  by  lateral  vortices.  1 — Total  induced  velocity. 


Table  l  The  harmonic  components  of  the  induced  Telocity  vf 
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Table  2  The  harmonic  components  of  the  induced  Telocity 


r  I 

t  X  l#'1  I  fjji,  *  I#"*  vsi/XlO**  j  o2W*10'* 


■  WIHMBM 


If*  P*V 


1 


REFERENCES 


[1]  [U.S.S.R]  Mili/  et  el.  Computation  and  Design  of  Helicopter 
(vol.  1/  Aerodynamics),  National  Defense  Industry  Publishing  Company, 
(1977),  p.  178. 

[2]  Wang  shicun,  The  Induced  Velocityof  Lift  Rotor,  Journal  of  North¬ 
western  Industrial  University,  (1963)  January. 

[3]  Wang  Zhuxi,  Guo  Duenren,  Introduction  of  Special  Functions,  (1965), 
p.  152. 

C 4  3  James  Scheitnaa,  A  Tabatatioa  of  Helicopter  Rotor-Blade  Differential  Pressures.  Stresses  and  Mo¬ 
tions  as  Measured  is  Flight.  NASA  TM  X-952,  Langley  Research  Center.  Langley  Station.  Ha¬ 
mpton.  Va..  December  tl.  (IMS). 


1 


Vi 


S 


SOME  PROBLEMS  IN  THE  FINITE-DIFFERENCE 
COMPUTATION  OF  THREE-DIMENSIONAL 
TRANSONIC  FLOWS 

Chen  Tiembt 

( Beijing  Institute  of  Aerodynamics ) 

Abstract 

The  mixed  finite-difference  relaxation  iteration  method  is  applied  to 
calculation  of  the  wing-body  combination  with  rectangular  wings  based  on 
the  three-dimensional  transonic  small-disturbance  potential  equation  in 
the  cylinder  coordinates.  Meanwhile*  the  influences  of  different  computa¬ 
tion  regions  and  relaxation  parameters  (on  subsonic  points)  on  the  ranges 
of  the  computed  Mach  numbers*  angles  of  attack  and  calculated  results 
are  studied. 

It  is  shown  that  extending  the  computation  region  brings  about  inc¬ 
reasing  the  convergence  range  of  the  small-disturbance  equation,  and  the 
convergence  rate  and  computation  accuracy  would  be  are  enhanced  if  the 
subsonic  relaxation  parameter  is  taken  to  be  1.9. 

Also*  it  is  demonstrated  that  the  calculated  results  differ  when  the 
finite  or  infinite  region  is  taken  as  the  computation  region  of  the  r-direc- 
tion.  Particularly  the  distribution  of  the  pressure  cooeficients  near  the 
wingtip  varies  obviously  as  the  incident  Mach  number  becomes  larger. 


SOME  PROBLEMS  IN  THE  FINITE-DIFFERENCE  COMPUTATION  OF  THREE- 
DIMENSIONAL  TRANSONIC  FLOWS 

Chen  Tiemin 

Beijing  Institute  of  Aerodynamics 
Submitted  25  January  1985 
I.  Introduction  of  Basic  Method 

1.  Three-Dimensional  Transonic  Small-Disturbance  Equation 

The  three-dimensional  transonic  small-disturbance  equation  in 
cylindrical  coordinates  is. 

AV„+<p„+-±-vr  +  -*.  qp«-  0  (D 

where  <p  is  disturbance  potential;  A  *P*—(  Y  +  1  )A/2.<P.»  3*— l  —  y 
is  specific  heat;  M\m  is  the  incident  Mach  number.  If  r-direction 
is  selected  as  infinite  region/  through  coordinate  conversion 
* —br/(l +ar),  the  infinite  region  of  r  is  converted  to  the  finite 
computation  region  /J  .  Here  a  and  b  are  constants  to  be  determined. 
Now  Equation  (1)  can  be  rewritten  as 


+_4i£L»„]+liL-I!)1*..,  (2) 

Then  Equations  (1)  and  (2)  correspond  to  the  finite  and  infinite 
regions  (select  proper  a  and  b  values)  of  r  respectively. 


2.  Boundary  Conditions  (See  Fig.  1) 


Fig.  1.  Boundary  surfaces  and  diagrammatic  drawing  of  model 
Key:  (l)Vertical  plane  of  symmetry;  (2)  Upstream  plane;  (3)  Outer 

cylindrical  surface;  (4)  Horizontal  plane;  (5)  Downstream  plane; 
(6)  Tail-trace  plane. 


On  the  surface  of  the  body:  Assume  the  surface  equation  of  the 
body  symmetrical  to  the  axes  is  r=R(x),  then  the  boundary  condition 
of  the  tangential  flow  is 

IVTL*r<P,'^“*(  *)- a  sin  0  3  (3) 
where  o<  is  the  incident  attack  angle. 


On  the  surface  of  the  wing:  Assume  the  equations  of  upper  and 

lower  wing  surfaces  are  Z —/*(*.  y)  ,  then  the  flow  along  the  upper 

and  lower  wing  surfaces  should  satisfy 

lim  (<P»)- r  C/i(  *,  y)-aj  (4) 

•  -«i 

Here  /«(*♦  >)“*/»(*,  0£  represents  the  upper  and  lower  wing  sur¬ 

faces  . 

On  the  tail-trace  surface:  It  is  known  from  the  Kutta  condition 

that 

v(x,  y,  oo-tpc*,  o.)-r(y)  (5) 

There  are  two  conditions  on  the  upstream  plane  and  the  outer 
cylindrical  surface: 

(1) If  Equation  (1)  is  used,  it  should  satisfy,  from  the  solutions 
of  small-disturbance  wing  lift  problems, 

where  a  —  y/'xi + ,  t  is  half  of  the  span  length. 

(2)  If  Equation  (2)  is  used,  the  disturbance  potential  is  0, 
i  .  e . , 


There  are  also  two  conditions  on  the  downstream  plane: 


(1)  If  Equation  (1)  is  used 


should  satisfy 
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where 


P-K  r1— 2rycos  8  +>*  . 


(2)  If  Equation  (2)  is  used,  the  disturbance  velocity 

<P.-0  (< 


3.  Finite-Difference  Formula,  Finite-Difference  Equation  and  Their 
Solutions 


This  paper  uses  the  variable  mixed  finite-difference  formula 
of  Murman-Cole,  i.e.,  C?  at  subsonic  point  uses  central  finite- 

J  XX 

difference;  at  supersonic  point  uses  incident  finite-difference; 
r  (or  )  and  0-direction  both  use  central  finite-difference.  The 
finite-difference  formula  at  special  point  (with  boundary  conditions 
substituted ) : 


(1)  According  to  boundary  condition  (3),  it  can  be  obtained  on 


the  surface  of  the  body  that: 
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(2)  On  the  upper  and  lower  wing  surfaces 
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(3)  On  the  tail-trace  plane,  the  pressure  continuum  and  Kutta 
Condition  are  applied 

]  ( 1 2 1 

Using  the  above  finite-difference  formulae,  and  on  the  plane 
of  i=constant  and  the  straight  live  of  k=constant,  the  algebraic 
simultaneous  equations  with  O  .  .  ,  .  ,  (Q .  .  .  and  CD.  .  .  .  as  un- 

*  lfj— 1/K  /  1  /  J  /  K  /  1  ,  J  +  1  .  K 

knowns  are  established.  These  equations  are  three-diagonal  and  can 

A 

be  solved  by  the  chase  method.  Assume  the  solution  is  then 

the  solution  after  relaxation  is 

«pr‘-<D$7‘  +  (l  (13) 

where  to  is  the  relaxation  parameter.  The  subsonic  relaxation  param¬ 
eter  eot  in  the  computation  is  selected  as  1.9  or  0.9;  the  supersonic 
relaxation  parameter^is  selected  as  0.9.  The  iteration  method  is 
conducted  to  solve  for  the  entire  flow  field.  The  iteration  stops 
when  the  entire  flow  field  satisfies 

d4) 

where  £  is  the  convergence  range  and  is  generally  selected  as 
£  =0. lxlO*4. 

The  computation  formula  for  pressure  coefficient  on  wing  is 

C,--2V.  (15) 

On  the  surface  of  the  body,  the  second-order  term  of  the  horizontal 
disturbance  velocity  is  reserved 

c,~, r*l)  (16) 
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II.  Example 


s 


m 


1.  Wing-Body  Combination  Equation 


(1)  The  wing  is  a  rectangular  wing  with  double-arc  airfoil 

Z“'T"(~7  *)  (17) 

The  airfoil  relative  thickness  ratio  ^=0.06;  chord  length  b=l  is 
a  reference  length;  the  aspect  ratio  ^_=L/b=5.4;  L  is  the  span  length 
including  the  body. 


2.  The  Body  is  the  Model  "AGARD-B"  Body,  the  Equation  for  the  Contour 
of  the  Head  is 

/?(*) - 1-  (5D+*)[l-  (5Z?+*)*  +  ~4-(5D+*)*] 


(  18) 


(-SD<*<2D) 

where  D  is  the  diameter  of  the  cylindrical  portion  of  the  body. 
Assume  D=l,  when  x  }  2 ,  then  R=0.5  and  extends  into  the  downstream. 


The  origin  of  coordinates  is  selected  at  the  intersection  of 
the  line  connecting  the  center  points  of  the  chord  and  the  body  axes, 
The  finite  region  of  r  is  marked  as  I  and  the  infinite  region  as 
II.  The  sampling  points  in  the  x,  r  (or  ^  )  and  9-direction  are 
all  distributed  at  nonequal  distance. 


III.  computation  Results 


1.  The  maximum  Mach  number  and  attack  angle  in  the  computation 


region  I  are  M  =0.926  and  CX  =2  respectively;  in  the  computation 

3  W 


region  II  they  are  Mapo=0.98  and  $  =4  respectively.  This  indicates 
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that  the  coordinate  conversion  ^=br/(l+ar)  in  the  r-direction  can 
increase  the  range  of  the  calculated  Mach  number  and  attack  angle 
of  the  small-disturbance  equation. 


2.  Figure  2  gives  the  comparison  curves  between  the  computed 

results  and  the  two-variable  test  results  at  M  .*.=0.965,  CK  =  2°  and 

a  cai 

y=0.5  (y  is  the  spanwise  relative  location  on  the  wing  cross-section), 
This  figure  shows  that  the  conmputed  results  in  region  II  are  reason¬ 
able  . 


( 1 ) 


"f«(  2) 


- (3) 

•  •••  =*Stt<M(4) 

Fig.  2.  Comparison  between  computed  results  and  test  results 
Key:  (1)  upper  surface;  (2)  lower  surface;  (3)  computed  results; 

(4)  two-dimensional  test  results. 


3.  Figure  3  shows  that  there  is  difference  between  the  computed 

results  of  region  I  and  II.  At  £  =0.1x10”^,  when  the  M  is  larger, 

there  is  a  difference  in  the  chordwise  distribution  of  wing  surface 

pressure  coefficients  and  it  becomes  more  obvious  near  the  wingtip. 

This  is  because  the  outer  boundary  of  region  I  is  like  a  wall,  causing 

the  local  Mach  number  to  increase  and  the  absolute  value  of  C  on 

P 

the  upper  surface  to  be  larger  and  that  for  the  lower  surface  to 
be  smaller. 


;i!< 
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van,' 


/Wo„-0.IS8,«*i  ,y»9.J 


Fig.  3.  Influence  of  computation  regions  on  computed  results 
Key:  (1)  upper  wing  surface;  (2)  lower  wing  surface. 


4.  Figure  4  gives  the  curves  of  influence  of  subsonic  relaxation 
parameters  on  computed  results.  Under  the  same  £,  the  computed 
results  at  &>^=1. 9  is  more  accurate  than  that  at  CO ^=0.9;  its  conver¬ 
gence  rate  is  faster.  There  are  similar  phenomena  in  both  computation 
region  I  and  II.  This  indicates  that,  under  guaranteed  convergence 
condition,  the  larger  CO  ^  is  selected  (>  2),  the  faster  the  convergence 
rate  of  computation  and  the  higher  the  accuracy  of  data. 


/  - J. »,«•!»-«,  A»J» 

0.4  '  — •,»o.*,«»10*\A""'n4 

,  •  ' 

Fig.  4.  Influence  of  relaxation  parameters  on  computed  results 
Key:  (1)  upper  wing  surface;  (2)  lower  wing  surface. 
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Abstract 


A'  three-dimensional  computation  of  a  jet  in  a  crossflow  has  been  per¬ 
formed  on  a  microcomputer.  The  SIMPLE  method  is  adopted  with  some 
modifications.  A  special  computer  code  is  developed  in  conformity  with  the 
limited  memories  of  the  microcomputer.  Velocity  and  presure  distributions 
are  '£fven  and  compared  with  the  available  experimental  data  satisfacto¬ 
rily. 
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Pan  Huachen  and  Zhang  Shiying 
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A  vertical  jet  shot  out  of  a  flat  plate  will  interfere  with  a 
low  speed  incident  flow  parallel  to  the  flat  plate  (crossflow),  caus¬ 
ing  complex  flow  phenomena  to  occur.  The  jet  is  deflected  by  the 
crossflow,  and  its  velocity  is  rapidly  reduced;  the  incident  flow 
is  blocked  and  drawn  by  the  jet  (Fig.  1).  The  jet  and  the  incident 
flow  interreact  with  each  other  causing  two  symmetrical  vortex  cores 
to  develop  under  the  jet. 


Fig.  1.  Flow  field  of  a  vertical  jet  interferring  with  a  crossflow 

and  its  computation  domain 

Key:  (1)  crossflow;  (2)  vortex. 

The  vortex  caused  by  this  interference  develops  along  the  direction 

of  the  wall  surface  and  it  can  be  used  as  a  means  to  control  the 

attach  surface  layer^1^.  It  can  also  be  applied  to  heat  transfer 

engineering  and  other  areas.  There  have  been  numerous  experimental 
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studies  on  this  kind  of  flow  both  in  our  country  and  abroad1- 
Due  to  the  difficulties  in  measurement,  it  is  generally  difficult 
to  obtain  detailed  data  of  the  entire  flow  field,  especially  near 
the  nozzle.  It  is  impossible  to  conduct  computation  using  two-dimen¬ 
sional  or  nonviscous  method  because  of  the  complexity  of  the  flow; 
the  full  three-dimensional  viscous  numerical  computation  method  must 
be  used.  This  kind  of  flow  possesses  strong  three-dimensional  mixing 
effects  and  needs  to  use  appropriate  turbulence  model.  According 
to  the  principles  in  References  [6,  7],  we  developed  a  computer  program 
and  performed  computation. 


I.  Basic  Equations 


The  basic  equation  of  a  three-dimensional,  steady,  incompressible 
time-averaged  Reynold  number  turbulent  flow  in  the  rectangular  coordi 
nates  is:  continuity  equation 


_£Hi. 

ax, 


(1) 


The  momentum  equation  (with  the  assumption  of  vortex  viscosity) 

'  ax,  ax,  ax,  1/ \  ax,  ax,  J J 
where  the  effective  viscosity  jbL  is  defined  as 

(3) 

where  pi  is  the  laminar  viscosity.  JJ.  is  the  turbulent  viscosity 
and  is  determined  by  the  following  turbulence  model 

v,-C.PK  7<  (4) 

where  K  is  the  turbulence  kinetic  energy.  €■  is  the  dissipation  rate 
of  K  and  is  determined  by  the  following  equations 

p,'i£— £r(-Sr  i£)+c-p«  (5) 


p“'7$-m-£r(-Sr  it)+(CiG-c'P()ir 


(6) 


where  G  is  the  generation  rate  of  the  turbulence  kinetic  energy 


The  coefficients  used  in  the  equations  are  adopted  from  the  data 

recommended  by  Reference  [8]: 

C.-0.09,  C.-1.44,  Cj-1.92,  a»-l,  o.-l.J  (8) 

Since  the  above  turbulence  model  is  used,  p*  in  Equation  (2)  is 

j>*-  p  +— |-P K 

where  2/3  K  can  be  regarded  as  turbulent  pressure  term. 


II.  Solution  Method 

The  method  used  by  this  paper  is  based  primarily  upon  the  SIMPLE 
r  7  I 

method  .  Said  method  adopts  the  iteration  solution  method.  The 
three  velocity  components  u,  v  and  w  used  to  solve  the  x-direction 
momentum  equation  are  all  old  values;  u  used  to  solve  the  y-direction 
momentum  equation  is  a  new  value;  and  u  and  v  used  to  solve  the  z- 
direction  momentum  equation  are  new  values.  It  was  found  during 
actual  computation  that  this  kind  of  nonconforming  value  character 
was  an  important  factor  which  caused  divergence  in  solving  the  equa¬ 
tion.  Therefore,  when  we  were  developing  the  program  the  newly  obtain¬ 
ed  u  and  v  values  were  temporarily  stored  and  old  u,  v  and  w  values 
were  used  to  perform  iteration  for  all  three  momentum  equations, 
thereby  guaranteeing  the  conforming  value  character. 

In  order  to  increase  efficiency  and  lower  cost,  the  existing 
INTEL-86-330  microcomputer  was  adopted  to  perform  computation.  Said 
program  requires  more  memory.  The  example  in  this  paper  requires 
more  than  100  giga  bytes  in  memory  at  least,  and  the  available  memory 
space  of  the  said  microcomputer  is  only  a  little  more  than  100  K. 
Therefore,  during  the  development  of  the  program  the  technique  of 
exchanging  with  secondary  storage  was  adopted.  Every  three-dimensional 
group  was  divided  into  a  dozen  or  so  two-dimensional  group  and  stored 
in  hard  disk.  During  computation  they  were  called  in  separately 
by  a  subroutine  to  take  part  in  the  iteration  and  then  restored. 

This  greatly  saved  memory  space.  Yet  as  a  price,  the  computation 
speed  dropped  about  20  fold. 


III.  Example 


The  experimental  conditions  in  Reference  [3]  for  the  incident 
flow  and  geometric  data  are  used.  The  incident  flow  velocity  is 
39  m/sec  and  jet  velocity  is  156  m/sec  with  a  velocity  ratio  of  4. 

The  symmetry  of  flow  field  only  considers  half  of  the  flow  field 
cut  by  the  symmetric  plane,  as  shown  in  Fig.  1.  The  computation 
domain  is  a  box-shaped  region  with  its  bottom  connected  to  the  flat 
plate.  One  of  the  sides  is  the  symmetric  plane  and  the  other  four 
sides  are  the  upstream,  downstream  boundaries  and  infinite  boudnaries, 
respectively.  The  grid  is  in  the  form  of  rectangular  coordinates: 
x-direction  (flow  direction)  has  20  rows  of  grid,  z-direction  (per¬ 
pendicular  to  the  plate)  has  16  rows  and  y-direction  has  12  rows 
for  a  total  of  3840  grid  points  (measurement  points).  In  order  to 
properly  simulate  the  details  of  flow  field  near  the  jet  nozzle, 
a  variable  spacing  grid  is  used  to  give  dense  grid  near  the  jet  nozzle 
and  sparse  grid  for  further  flow  field.  Using  the  nozzle  diameter 
D  as  the  reference,  the  upstream  boundary  is  3.5  D  from  the  center 
of  the  nozzle  and  the  downstream  boundary  is  8  D  from  the  center 
of  the  nozzle.  The  z-direction  infinite  boundary  is  about  12  D  from 
the  plate  and  the  side  infinite  boundary  is  about  5  D  from  the  symmet¬ 
ric  plane.  The  method  in  Reference  [6]  is  adopted.  Six  grid  points 
are  used  to  simulate  a  semi-circle  jet  nozzle.  Meanwhile  the  total 
area  corresponding  to  the  six  grid  points  is  made  to  be  in  agreement 
with  the  area  of  the  semi-circle  jet  nozzle. 


The  upstream  boundary  defines  the  incident  flow  velocity.  The 
normal  velocities  on  the  infinite  boundary,  the  symmetric  plane  bound¬ 
ary  and  the  flat  plate  are  all  assumed  to  be  zero.  The  jet  velocity 
is  defined  by  the  jet  nozzle.  The  gradient  of  the  tangential  veloci¬ 
ties  v  and  w  on  the  downstream  boundary  are  assumed  to  be  zero,  yet 
u  is  adjusted  according  to  the  total  flow  rate  on  the  basis  of  its 

value  on  the  previous  row.  The  wall  surface  is  processed  by  adopting 

r  Q  I 

the  wall  surface  function  method1  . 

The  standard  for  computation  convergence  is  set  at  when  both 
the  corresponding  remainders  of  the  momentum  and  continuity  equations 
reach  5%.  This  example  required  90  steps  and  spent  170  hours  on 
the  microcomputer. 

V.  Results  and  Their  Discussions 

In  order  to  verify  the  computation,  the  results  are  compared 

with  the  test  data  of  Reference  [3,  4,  5].  Figure  2  shows  the  velocity 

contours  on  the  symmetric  plane.  In  the  figure,  the  line  with  indexing 

points  is  the  jet  center  line  drawn  by  using  the  peaks  of  the  contours, 

and  it  can  be  seen  that  the  line  is  in  good  agreement  with  the  test 

* 

data.  See  Fig.  3  for  the  decay  of  jet  velocity.  In  the  figure, 

W  is  the  jet  velocity  at  the  jet  nozzle,  s  is  the  natural  coordi- 
nates  along  the  jet,  is  the  jet  velocity  along  s.  The  computed 
velocity  decays  faster  than  the  test  value,  but  the  difference  is 
not  distinctive  and  both  show  that  the  rate  of  velocity  decay  follows 
a  rule  of  first  small,  then  large,  and  again  from  large  back  to  small. 
Figure  4  is  the  projection  of  vortex  core  contour  on  the  symmetric 


plane.  The  vortex  core  trace  derived  from  the  computed  values  is 
more  complicated.  The  low  pressure  centers  basically  coincide  with 
the  vortex  core  contour,  otherwise  the  centripetal  force  that  makes 
the  air  streamlines  turn  vortically  will  not  be  generated.  The  curve 
in  Fig.  4  is  the  projection  of  the  low  pressure  centers  on  the  symmet 
ric  plane,  i.e.,  the  approximation  of  the  vortex  core  contour,  and 
matches  well  with  the  vortex  core  test  values.  Figure  5  is  the  com¬ 
parison  of  computed  values  (left  half)  and  test  values  (right  half) 
of  the  pressure  coefficient  distributions  on  the  flat  plate.  Judging 
from  the  trend  of  pressure  distributions,  the  computed  values  of 
pressure  rise  behind  the  jet  nozzle  on  the  flat  plate  are  faster 
than  the  test  values. 


Fig.  2.  Velocity  distributions  and  the  location  of  jet  center  line 
on  the  symmetric  plane 

Key:  (1)  test  values  of  Reference  [3];  (2)  test  values  of  Reference 
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Fig.  3.  Rule  of  velocity  decay  on  the  jet  center  line 
Key:  (1)  test  values  of  Reference  [4], 


Fig.  4.  Projection  of  the  location  of  vortex  core  on  the  symmetric 
plane 

Key:  (1)  test  values  of  Reference  [31. 


3 •  Pressure  coefficient  distribution  on  the  flat  plate 

(1)  Crossflow;  (2)  Computed  values;  (3)  test  values  of  Reference 
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The  three-dimensional  viscous  flow  computation  method  introduced 
by  this  paper  can  compute  detailed  data  of  the  interference  flow 
field  of  a  vertical  jet  in  an  incident  flow.  The  comparison  with 
test  results  shows  that  the  computed  results  are  more  accurate,  espe¬ 
cially  the  locations  of  jet  contour,  vortex  core  contour,  etc.  The 
fact  that  the  computed  jet  decay  is  a  little  faster  and  the  computed 
pressure  rise  behind  the  jet  nozzle  is  faster  indicates  that  the 
numerical  dissipation  might  be  larger,  but  not  enough  to  cause  signif¬ 
icant  influence  on  the  prediction  of  the  entire  flow  field.  Therefore 
the  said  numerical  method  can  be  an  effective  tool  for  analyzing 
flow  field. 
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